home *** CD-ROM | disk | FTP | other *** search
/ MacFormat 1994 August / August CD.bin / Shareware / Education / numericalmethods Folder / chap_6 / a6_3.m < prev    next >
Encoding:
Text File  |  1994-06-05  |  3.6 KB  |  154 lines  |  [MATS/MATL]
  1. echo off;
  2. % NUMERICAL METHODS: MATLAB Programs, (c) John H. Mathews 1994
  3. % To accompany the text:
  4. % NUMERICAL METHODS for Mathematics, Science and Engineering, 2nd Ed, 1992
  5. % Prentice Hall, Englewood Cliffs, New Jersey, 07632, U.S.A.
  6. % This free software is complements of the author.
  7.  
  8. % Algorithm 6.3 (Differentiation Based on N+1 Nodes).
  9. % Section    6.2,    Numerical Differentiation Formulas, Page 342
  10. echo on; clc; format long; hold off; clear
  11.  
  12. % This program differentiates a Newton polynomial approximation.
  13. % n+1  are points needed to construct  Pn(x).
  14.  
  15. % The abscissas are stored in  X.
  16. % The ordinates are stored in  Y.
  17. % The points are counted k = 1,2,...,n+1.
  18.  
  19. % The example investigates f(x) = cos(x) using 5 nodes.
  20.  
  21. % The derivative is approximated at X(1).
  22.  
  23. %    The function f(x) is stored in the M-file  f.m.
  24. % function z = f(x)
  25. % z = cos(x);
  26.  
  27. delete f.m
  28. diary f.m; disp('function z = f(x)');...
  29.            disp('z = cos(x);');...
  30. diary off;
  31.  
  32. f(0); % Remark. f.m and diffnew.m are used for Algorithm 6.3
  33.  
  34. pause % Press any key to continue.
  35.  
  36. clc;
  37. %    Place abscissa for differentiation in x0.
  38.  
  39. % Place the endpoints interval [a0,b0] about x0 in a0 and b0.
  40.  
  41. x0 = 0.8;
  42.  
  43. a0  = 0;
  44.  
  45. b0  = pi/2;
  46.  
  47. pause    % Press any key to graph the function.
  48.  
  49. clc; clg;
  50. y0 = f(x0);
  51. hs = (b0-a0)/100;
  52. Xs = a0:hs:b0;
  53. Ys = f(Xs);
  54. a = 0;
  55. b = 1.6;
  56. c = 0;
  57. d = 1.1;
  58. axis([a b c d]);...
  59. plot(x0,y0,'o',Xs,Ys);...
  60. hold on;...
  61. plot([a b],[0 0],'b',[0 0],[c d],'b');...
  62. xlabel('x');...
  63. ylabel('y');...
  64. title('y = f(x) and the given point');...
  65. grid;...
  66. axis;...
  67. hold off;...
  68. shg; pause    % Press any key to continue.
  69.  
  70. clc;
  71. %    Place abscissa for differentiation in  x0.
  72.  
  73. % Place the number of points in  n.
  74.  
  75. % Place the step size in  h.
  76.  
  77. x0 = 0.8;
  78. n = 5;
  79. h = 0.01;
  80. X = x0:h:x0+(n-1)*h;
  81. Y = f(X);
  82.  
  83. [A,df] = diffnew(X,Y);
  84.  
  85. pause % Press any key to continue.
  86.  
  87. format short;
  88. Mx0 = 'Approximating the derivative by differentiating a Newton polynomial.';
  89. Mx1 = 'First construct the Newton approximation polynomial.';
  90. Mx2 = 'The abscissas are:';
  91. Mx3 = 'The ordinates are:';
  92. Mx4 ='The coefficients for the Newton polynomial are:';
  93. Mx5 = 'The centers for the Newton polynomial are:';
  94. clc,echo off,diary output,...
  95. disp(''),disp(Mx0),disp(''),disp(Mx1),disp(''),disp(Mx2),disp(X),...
  96. disp(Mx3),disp(Y),disp(''),disp(Mx4),format short,disp(A),...
  97. disp(Mx5),disp(X(1:n-1)),diary off, echo on,format long
  98.  
  99. pause    % Press any key to see the approximate derivative.
  100.  
  101. clc;
  102. X1 = [a b];
  103. C = [df f(x0)];
  104. Y1 = polyval(C,X1-x0);
  105. axis([a b c d]);...
  106. plot(x0,y0,'or',Xs,Ys,'-g',X1,Y1,'-r');...
  107. hold on;...
  108. plot([a b],[0 0],'b',[0 0],[c d],'b');...
  109. xlabel('x');...
  110. ylabel('y');...
  111. title('The tangent line has slope m = f`(x0).');...
  112. grid;...
  113. axis;...
  114. hold off;...
  115. shg; pause    % Press any key to continue.
  116.  
  117. Mx1 = 'Numerical approximation for the derivative.';
  118. Mx2 =  ['f`(',num2str(X(1)),')'];
  119. Mx3 =  ['P`(',num2str(X(1)),')'];
  120. Mx4 = [Mx2,' ~ ',Mx3,' = '];
  121. clc,disp(''),disp(Mx1),...
  122. disp(''),disp(Mx4),disp(df)
  123.  
  124. pause    % Press any key to zoom in.
  125.  
  126. a = min(X)-3*h;
  127. b = max(X)+3*h;
  128. c = min(Y)-3*h;
  129. d = max(Y)+3*h;
  130. X1 = [a b];
  131. C = [df f(x0)];
  132. Y1 = polyval(C,X1-x0);
  133. axis([a b c d]);...
  134. X0 = X(2:n);...
  135. Y0 = Y(2:n);...
  136. plot(x0,y0,'or',X0,Y0,'og',Xs,Ys,'-g',X1,Y1,'-r');...
  137. hold on;...
  138. plot([a b],[0 0],'b',[0 0],[c d],'b');...
  139. xlabel('x');...
  140. ylabel('y');...
  141. title('The tangent line has slope m = f`(x0).');...
  142. grid;...
  143. axis;...
  144. hold off;...
  145. shg; pause    % Press any key to continue.
  146.  
  147. Mx1 = 'Numerical approximation for the derivative.';
  148. Mx2 =  ['f`(',num2str(X(1)),')'];
  149. Mx3 =  ['P`(',num2str(X(1)),')'];
  150. Mx4 = [Mx2,' ~ ',Mx3,' = '];
  151. clc,echo off,diary output,...
  152. disp(''),disp(Mx1),...
  153. disp(''),disp(Mx4),disp(df),diary off,echo on
  154.